function myDist_new = update_dist_int_edu_ss(myDist, Pr, Ss, rs)
    global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha xi xi1 ub lb pi_e...
        w theta Int_pi_w sigmae elasts log_e_mu log_e_sigma eta phi sigma gamma rts w_mult wspread unemp
    
    myDist_new = zeros(size(myDist));
    ln_e_hat = -1/2*log_e_sigma^2;
    [~, F_es, wprimeS, ~] = choiceProb(Ss, rs);

    w_bar_p = sum(myDist .* wprimeS, 1:3);
 
    w2 = w_bar_p * wspread;
    w_max = max(w2);
    w_min = min(w2);
 

    for c = 1:length(Ss)
        for j = 1:N_a
            for i = 1:N_w
                wage = w(i);

                t1 = unemp + (b   + squeeze(F_es(i,j,c))')  .* f_w(wage);
                midpts = log(w2) - log(t1);
                pts = ([midpts(1,:);(midpts(2:end,:) + midpts(1:end-1,:))/2; Inf] - ln_e_hat)/log_e_sigma;
                mean_prob = diff(normcdf(pts));
                tmp = bsxfun(@times, mean_prob, Pr(j,:));
                pr1 = myDist(i, j, c) .*  tmp ;
                myDist_new(:, :, c) = myDist_new(:, :, c) + pr1;
                
            end
        end
    end
       
    w = w2;
end